A Multilevel Blocking Approach to the Sign Problem in Real-Time Quantum Monte 

Carlo Simulations 
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We propose a novel approach toward the general solution 
of the sign problem in real-time path-integral simulations. 
Using a recursive multilevel blocking strategy, this method 
circumvents the sign problem by synthesizing the phase can- 
cellations arising from interfering quantum paths on multiple 
timescales. A number of numerical examples with one or a 
few explicit degrees of freedom illustrate the usefulness of the 
method. 



Path integrals [[!) provide an elegant alternative to the 
operator formulation of quantum mechanics. Because 
they are easily adapted to many-body systems, quantum 
Monte Carlo (QMC) simulations based on path integrals 
can potentially yield exact results for the dynamics of 
condensed phase quantum systems. A number of early 
attempts to use QMC simulations for real-time dynamics 
demonstrated their potential, but these studies also 
uncovered the ubiquitous "dynamical sign problem" - 
interference among quantum paths leads to large statisti- 
cal noise that increases linearly with the number of pos- 
sible paths, which in turn grows exponentially with the 
timescale of the problem. Consequently, real-time QMC 
simulations have been limited to problems of very short 
timescales. Several attempts to extend the timescale of 
real-time QMC simulations have appeared |5j-[|, all of 
which rely on a common idea — by using various filter- 
ing schemes, the phase cancellations can be numerically 
stabilized by damping out the non-stationary regions in 
configuration space. Such filtering methods were able to 
extend the timescale somewhat, but they were generally 
unable to reach timescales of interest in typical chemical 
systems. 

Here we propose a novel approach based on a blocking 
strategy which may provide a general solution to the dy- 
namical sign problem for very long times. The blocking 
strategy asserts that by sampling groups of paths called 
blocks, the sign problem is always reduced compared to 
sampling single paths ||. Though this blocking strat- 
egy seems simple, its practical implementation is cum- 
bersome, especially when going out to long times. Be- 
cause the number of paths grows exponentially with the 
timescale, the number of blocks also grows immensely. To 
cure this problem, we first define elementary blocks and 
then group them together into larger blocks. Blocks of 



different sizes are introduced on several timescales called 
levels. After taking care of the sign cancellations within 
all blocks on a finer level, the entire sign problem can then 
be transferred to the next coarser level. By recursively 
proceeding from the bottom level (shortest timescale) up 
to the top (longest timescale), the troublesome numeri- 
cal instabilities associated with the sign problem can be 
systematically avoided. In a slightly different form, this 
multilevel blocking (MLB) algorithm has recently been 
applied to treat the "fermion sign problem" in many- 
fermion imaginary-time simulations [ ^] . As most techni- 
cal details of the MLB scheme can be found in , we only 
give a brief description of the algorithm below, stressing 
the main ideas and its differences from the fermion for- 
mulation. 

For concreteness, we describe the MLB method for the 
calculation of an equilibrium time-correlation function, 
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With minor modifications, the MLB method also ap- 
plies to other dynamical properties like the ther- 
mally symmetrized correlation function j§], C s (t) = 

Z -l Tl { e -(ph/2+U)H/h Ae -((3/K/2-U)H/h B ^ ^ % bdng 

the partition function. In terms of path integrals, the 
traces in ([!]) involve two quantum paths, one propagated 
backward in time for the duration —t and the other prop- 
agated in complex time for the duration t — i/3h. Dis- 
cretizing each of the two paths into P slices, the entire 
cyclic path has a total of 2P slices. A slice on the first 
half of them has length —t/P, and on the second half 
(t — i/3h)/P. We require P — 2 L which defines the total 
number of levels L. Denoting the quantum numbers (e.g., 
spin or position variables) at slice j by rj, {n, • ■ • r 2 p} 
is a discrete representation of a path, and the correlation 
function (\n) reads 



J dn ■ ■ ■ dr 2 pB(r2p)A(rp)Y\ 2 j ^ 1 (r J ,r j+1 )o 



J dn ■ ■ ■ dr 2P nf£iOi, r j+ i) 
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where the level-0 bond (rj, r 3+ i)o is simply the short-time 
propagator between slices j and j + I, and r 2 p +1 = T\. 
A direct application of the QMC method would sample 
these paths using the modulus of the integrand in the 
denominator as the weight. 

We first assign all slices along the discretized path to 
different levels i = 0, . . . , L (see Figure IT]). Each slice 
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j = 1 , . . . , 2P belongs to a unique level £, such that 
j = (2k + 1)2 and k is a nonnegative integer. For 
instance, slices j = 1,3,5,... belong to level £ = 0, 
slices j = 2,6,10,... to t = 1, etc. The MLB algo- 
rithm starts by sampling only configurations which are 
allowed to vary on slices associated with the finest level 
1 = 0, using the weight V = |(ri, r 2 ) • • • {r 2 p, ri) Q \. 
The short-time level-0 bonds are then employed to syn- 
thesize longer-time level- 1 bonds that connect the even-j 
slices. Subsequently the level-1 bonds are used to syn- 
thesize level-2 bonds, and so on. In this way the MLB 
algorithm moves recursively from the finest level (£ = 0) 
up to increasingly coarser levels until £ = L, where the 
measurement is done using r 2 p and rp. 
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FIG. 1. Example of how slices (circles) are assigned to 
various levels for L — 2. The first path goes from j = 8 to 
4 and the second from j = 4 to 8. The finest level 1 = 
contains j = 1,3,5,7, level I = 1 contains j = 2,6, and the 
top level 1 = 2 contains j = 4, 8. Coarse bonds are indicated 
by dotted (level-1) and dashed (level-2) lines. 

Generating a MC trajectory containing K samples for 
each slice on level £ = and storing these samples, we 
compute the level-1 bonds according to 

(r ] ,r ]+2 ) 1 = ($[(rj,r i+ i)o(r i+1 ,rj +2 )o]>-Po ( 3 ) 
= K' 1 $ [( r J> r J+i)o( r i+i. r j+2)o] , 

where the summation extends over the K samples, 

and <&[z] = e taTS ( z ' denotes the phase. For a complete 
solution of the sign problem, the sample number K has 
to be sufficiently large . The bonds (||) contain crucial 
information about the sign cancellations on the previous 
level £ = 0. Their benefit becomes clear when rewriting 
the integrand of the denominator in (||) as 

Vq x (r 2 ,r 4 )i • • • {r 2 p-2,r2p)i(r2P,r 2 )i . 

Comparing this to (^), we notice that the entire sign 
problem has been transferred to the next coarser level. 

In the next step, the sampling is carried out on level 
£ = 1 in order to compute the next-level bonds, using the 
weight VqPi with V\ = |(r2,r 4 ) 1 • • ■ (r 2 p,r 2 )i|. Gener- 
ating a sequence of K samples for each slice on level 1=1, 
and storing these samples, we then calculate the level-2 
bonds in analogy with (H) , 



{rj,r j+i ) 2 = ($[(rj,r i+2 )i(r j+2! ?- j+ 4)i]) Wl , 

and iterate the process up to the top level by employing 
analogously defined level-£ bonds. The correlation func- 
tion (Q) can then be computed from 

{B(r 2P )A(rp)^[(rp,r 2 p) L (r 2Pl rp) L \) v 
{®[(rp,r 2 p) L (r 2 p,rp) L ]) v 

with the positive definite MC weight V = VqP\ ■ ■ -Vl- 
The denominator in (Q) gives the average phase and in- 
dicates to what extent the sign problem has been solved. 
Under the direct QMC method, the average phase decays 
exponentially with t and is typically close to zero. With 
the MLB algorithm, however, the average phase remains 
close to unity even for long times, with a CPU time re- 
quirement ~ t. The price to pay for the stability of the 
algorithm is the increased memory requirement ~ K 2 as- 
sociated with having to store the sampled configurations. 

Now we illustrate the practical usefulness of the MLB 
method by several numerical examples. In each of these 
examples, we compute a time-correlation function. The 
average phase is larger than 0.6 for all data sets shown 
below. The decay in the average phase with t is a result 
of the finiteness of K . Choosing a larger K allows for 
a larger average phase out to longer time at the cost of 
increased computer memory and CPU time. Each data 
point in even the most intensive calculation took no more 
than a few hours on an IBM RS 6000/590. 

A. Harmonic oscillator. For H = p 2 /2m + muj 2 x 2 /2, 
the real and imaginary parts of (x(0)x(t)) oscillate in 
time due to vibrational coherence. In dimensionless units 
m = lu = 1, the oscillation period is 2tt. Figure |^(a) 
shows MLB results for C(t) = Re (x(O)x(t)). With P = 
32 for the maximum time t = 26, K = 200 samples were 
used for sampling the coarser bonds. Within error bars, 
the data coincide with the exact result and the algorithm 
produces stable results free of the sign problem. Without 
MLB, the signal-to- noise ratio was practically zero for 
t > 2. 

B. Two-level system. For a symmetric two-state sys- 
tem, H = ~^Aa x , the dynamics is controlled by tun- 
neling. The spin correlation function (a z (0)a z (t)) ex- 
hibits oscillations indicative of quantum coherence. Fig- 
ure |(b) shows MLB results for C{t) = Re(a z (0)a z (t)}, 
Putting A = 1, the tunneling oscillations have a period 
of 27r. With P = 64 for the maximum time t = 64, only 
K = 100 samples were used for sampling the coarser 
bonds. The data agree well with the exact result. Again 
the simulation is stable and free of the sign problem. 
Without MBL, the simulation failed for t > 4. 
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FIG. 2. MLB results (closed circles) for various systems. 
Error bars indicate one standard deviation, (a) C(t) for a 
harmonic oscillator at /3h = 1. The exact result is indicated 
by the solid curve, (b) Same as (a) for a two-level system 
at (3h = 10. (c) Same as (a) for a double-well system at 
/3h = 1. This temperature corresponds to the classical barrier 
energy, (d) C' a (t) for a double- well system coupled to two 
oscillators at f3% = 1. For comparison, open diamonds are for 
the uncoupled (a — 0) system. Note that C' B (t) is similar but 
not identical to C{i) shown in (c) . Solid and dashed lines are 
guides to the eye only. 



At low temperatures, interwell transfer occurs through 
tunneling motions on top of intrawell vibrations. These 
two effects combine to produce nontrivial structures in 
the position correlation function. At high tempera- 
tures, interwell transfer can also occur by classical bar- 
rier crossings. Figure ||(c) shows MLB results for C'(t) = 
Re (x(O)x(f)). The slow oscillation corresponds to inter- 
well tunneling, with a period of approximately 16. The 
higher-frequency motions are characteristic of intrawell 
oscillations. In this simulation, K — 300 samples were 
used. The data reproduce the exact result well, capturing 
all the fine features of the oscillations. Again the calcu- 
lation is stable and free of the sign problem, whereas a 
direct simulation failed for t > 3. 

D. Multidimensional tunneling system. As a final ex- 
ample, we consider a problem with three degrees of free- 
dom, in which a particle in a double-well potential is bi- 
linear ly coupled to two harmonic oscillators. The quartic 
potential in the last example is used for the double-well, 
and the harmonic potential in the first example is used 
for both oscillators. The coupling constant between each 
oscillator and the tunneling particle is a = 1/2 in dimen- 
sionless units. For this example, we computed the cor- 
relation function C s (t) for the position operator of the 
tunneling particle. Direct application of MC sampling to 
C s (t) has generally been found unstable for t > /3h/2 [0. 
In contrast, employing only moderate values of K = 400 
to 900, the MLB calculations allowed us to go up to 
t = 10/371. (Notice that this K is no larger than three, 
i.e. the number of dimensions, times the K needed for 
one-dimensional systems.) Figure 0(d) shows MLB re- 
sults for C' s (t) = ReC s (t). For the coupled system, the 
position correlations have lost the coherent oscillations 
and instead decay monotonically with time. Coupling to 
the medium clearly damps the coherence and tends to 
localize the tunneling particle. 

The data presented here demonstrate that the MLB 
method holds substantial promise toward an exact and 
stable simulation method for real-time quantum dynam- 
ics computations of many-dimensional systems up to 
timescales of practical interest. Instead of the exponen- 
tially vanishing signal-to-noise ratio in a ordinary appli- 
cation of the Monte Carlo method to real-time path inte- 
gral problems, the MLB method has a CPU requirement 
that grows only linearly with time. Moreover, the data 
we have so far seem to suggest that the memory require- 
ment K also grows only linearly with the dimensionality 
of the system. 
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C. Double-well potential. Next, we examine a double- 



well system with the quartic potential V(x) 
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